LDR to HDR

clear all; close all;

% Read image file

ldrImg1 = imread("C:\image processing lab\LDRtoHDR\LDR_022.jpg");

hdrImg1 = hdrread("C:\image processing lab\LDRtoHDR\HDR_022.hdr");

hdrImg1 = tonemap(hdrImg1);

LDRtoHDR(ldrImg1, hdrImg1);

The Loss Parameter: 17.2992

% Read image file

ldrImg2 = imread("C:\image processing lab\LDRtoHDR\LDR_018.jpg");

hdrImg2 = hdrread("C:\image processing lab\LDRtoHDR\HDR_018.hdr");

hdrImg2 = tonemap(hdrImg2);

LDRtoHDR(ldrImg2, hdrImg2);

The Loss Parameter: 17.974

% Read image file

ldrImg3 = imread("C:\image processing lab\LDRtoHDR\LDR_010.jpg");

hdrImg3 = hdrread("C:\image processing lab\LDRtoHDR\HDR_010.hdr");

hdrImg3 = tonemap(hdrImg3);

LDRtoHDR(ldrImg3, hdrImg3);

The Loss Parameter: 18.1217

function LDRtoHDR(ldrImg, hdrImg)

STEP 1

Pre - Processing: LDR Image Segmentationּ

    % Convert image to HSV color space
    I_hsv = rgb2hsv(ldrImg);
    
    % Extract the brightness value
    V = I_hsv(:,:,3);
    
    % Create histogram
    [counts, ~] = imhist(V);
    
    % Calculate Otsu's threshold
    level = graythresh(V);
    otsu_threshold = round(level * 255); 
    
    % Initialization
    valley_set = [];
    
    % Algorithm begins
    for pixel_value = otsu_threshold:252
        
        % Find the valley points
        slopes = [counts(pixel_value) - counts(pixel_value - 1),
                  counts(pixel_value) - counts(pixel_value - 2),
                  counts(pixel_value) - counts(pixel_value - 3),
                  counts(pixel_value + 1) - counts(pixel_value),
                  counts(pixel_value + 2) - counts(pixel_value),
                  counts(pixel_value + 3) - counts(pixel_value)];
              
        if all(slopes >= 0) % Valley point condition
    
            % Classify into two classes based on the valley point
            class1 = counts(1:pixel_value);
            class2 = counts(pixel_value+1:end);
    
            % Calculate the sum of variance
            s = var(double(class1)) + var(double(class2));
    
            % Include s and valley points in the valley set
            valley_set = [valley_set; pixel_value, s];
        end
    end
    
    % Find the valley point with the largest s and assign it as threshold
    [~, idx] = max(valley_set(:,2));
    thr = valley_set(idx,1);
    
    % Create mask
    mask = max(0, V - thr/(1 - thr)).^2;
    
    % Normalize mask
    mask = mask / max(mask(:));
    
    % Multiply the mask by the original image to segment it
    ldrImg = im2double(ldrImg);  % Convert I to double
    ImgBright = ldrImg .* mask;
    ImgDark = ldrImg .* (1 - mask);
    
    % Show results
    figure("Position",[0 0 1920 1080]);
    sgtitle('Step 1: Pre-Processing');
    subplot(1,3,1); imshow(mask); title('Mask');
    subplot(1,3,2); imshow(ImgBright); title('Bright Image');
    subplot(1,3,3); imshow(ImgDark); title('Dark Image');

STEP 2

Create the Brightness and Darkness Blocks

    % Calculate average brightness for the bright and dark images
    avgBright = mean(ImgBright(:));
    avgDark = mean(ImgDark(:));
    
    % Inverse proportional gamma values
    gammaBright = 1.0 / (avgBright + 0.1);  % Adding 0.1 to prevent division by zero
    gammaDark = 1.0 / (avgDark + 0.1);  % Adding 0.1 to prevent division by zero
     
    
    % Apply gamma correction
    ImgBright_gamma = ImgBright .^ gammaBright;
    ImgDark_gamma = ImgDark .^ gammaDark;
    figure("Position",[0 0 1920 1080]);
    sgtitle('Step 2: Brightness & Darkness Blocks');
    subplot(1,3,1); imshow(ldrImg); title('Original');
    subplot(1,3,2); imshow(ImgBright_gamma); title('Brighten Bright Image');
    subplot(1,3,3); imshow(ImgDark_gamma); title('Darken Dark Image');

STEP 3

Combine the 2 Images and Display Them

    % Combine the images
    img_gamma = ImgBright_gamma + ImgDark_gamma;
    
    % Display the result
    figure("Position",[0 0 1920 1080]);
    subplot(1,2,1); imshow(ldrImg); 
    subplot(1,2,2); imshow(img_gamma); 
 
    figure("Position",[0 0 1920 1080]);
    sgtitle('The Result');
    subplot(2,2,1); imshow(ldrImg); title('Original Image')
    subplot(2,2,2); imshow(img_gamma); title('Result');
    subplot(2,2,3); imhist(ldrImg); title('Histogram of the original');
    subplot(2,2,4); imhist(img_gamma); title('Histogram of the result');
    figure("Position",[0 0 1920 1080]);
    sgtitle('Compare to the real HDR image');
    subplot(1,3,1); imshow(ldrImg); title('Original LDR Image');
    subplot(1,3,2); imshow(hdrImg); title('Original HDR Image')
    subplot(1,3,3); imshow(img_gamma); title('My Result');

Loss Function

    % Assuming you have h_hat (predicted HDR image) and h (ground truth HDR image)
    epsilon = 1e-6;  % Small constant to prevent log(0)
    beta = 0.3;      % Loss function weight
    h_hat = double(img_gamma);
    h = double(hdrImg);
    LRGB_R = abs(log(h_hat(:,:,1) + epsilon) - log(h(:,:,1) + epsilon));
    LRGB_G = abs(log(h_hat(:,:,2) + epsilon) - log(h(:,:,2) + epsilon));
    LRGB_B = abs(log(h_hat(:,:,3) + epsilon) - log(h(:,:,3) + epsilon));
    LRGB = LRGB_R + LRGB_G + LRGB_B;
    
    % Compute LLuminance
    Y_hat = 0.2126 * h_hat(:,:,1) + 0.7152 * h_hat(:,:,2) + 0.0722 * h_hat(:,:,3);
    Y = 0.2126 * h(:,:,1) + 0.7152 * h(:,:,2) + 0.0722 * h(:,:,3);
    LLuminance = abs(log(Y_hat + epsilon) - log(Y + epsilon));
    
    % Compute total loss
    L = LRGB + beta * LLuminance;
    
    % Display loss
    disp(['The Loss Parameter: ', num2str(mean(L(:)))]);  % The mean function is used to compute the average loss per pixel
end